library cml;
#include cml.ext;
CMLset;

__title=" ** A Beta Model, Analytical Solution **";
dta="/home/gov/faculty/ppaolino/parr/gaydata1";
let dep=pvt;
let ind=task wealth rights elite col; @  wi; la de; @
let sel=0;

vars=dep|ind; 

dataset=listw(dta,vars,sel);

stval=    2.3243|
          0.0196|
         -2.2845|
	 -0.3056|
	  0.0123|
	  0.0123|
	 10.2340|
	 -0.2995|
	 -5.5561|
	 -0.9105|
	  0.0123|
	  0.0123;
      
proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

proc gradproc(theta,ds);
    local y,x,z,be,h,xb,p,q,per1,gr1,gr2,grad;
		     
    y=(ds[.,1]+1)*100/10101;	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    per1=p+q;
    gr1= psi(p+q).*x.*p-psi(p).*x.*p+x.*p.*ln(y);
    gr2= psi(p+q).*x.*q-psi(q).*x.*q+x.*q.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);
    endp;

proc loglik(theta,ds);
    local y,x,z,be,h,xb,p,q,l1,l2,l3,l4,l5,llik;
		     
    y=(ds[.,1]+1)*100/10101;
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    retp(llik);endp;

_cml_GradCheck=1;
_cml_GradProc=&gradproc; 
_cml_Options = { none }; 
_cml_DirTol=0.00001; 

    if ind==0;
      _cml_ParNames="const"|"const";
    else;
      _cml_ParNames="const"|ind|"const"|ind;
    endif;

t=dta $+ dep $+ ".std.out"; 
output file=^t on; 
    print "dependent variable=" ;; $dep; 
    {b,logl,g,vc,ret}=CMLprt(CML(dataset,0,&loglik,stval));
@ cond(_cml_FinalHess); @

/* set ind and indV matrices */
y=(dataset[.,1]+1)*100/10101;	     
one=ones(rows(dataset),1);

if ind==0;
  x=one;
  else;
  x=ones(rows(dataset),1)~dataset[.,2:rows(ind)+1]; 
endif;

/* get predicted values of y */
  b1=b[1:cols(x),1];
  h1=trimr(b,rows(b1),0);
  p1=exp(x*b1);
  q1=exp(x*h1);
  eyb=p1./(p1+q1);
  
/* set descriptives for x and z */
meanx=meanc(x);  @ these are k1X1 vectors @
stdx=stdc(x);	
minx=minc(x);	
maxx=maxc(x);	

/* create switching parameter */

if ind/=0;
sw=zeros(1,cols(x)-1)|eye(cols(x)-1);
xl=meanx-stdx.*sw;           @ these are kXk-1 vectors @
xh=meanx+stdx.*sw;       
xmin=meanx.*(1-sw)+minx.*sw;  
xmax=meanx.*(1-sw)+maxx.*sw;
endif;

al1=exp(xl'*b1);                          @ nxk-1 matrix @
bl1=exp(xl'*h1);
ah1=exp(xh'*b1);
bh1=exp(xh'*h1);

meanl1=al1./(al1+bl1);                    @ nxk-1 matrix @
meanh1=ah1./(ah1+bh1);

meanh1-meanl1;

/* generate the random parameter vector */
     betasim=(rndmn(b,vc,1000));  @ nxk matrix, where n=# of draws @
     be=betasim[.,1:cols(x)];                @ nxk matrix @
     h=betasim[.,cols(x)+1:cols(betasim)];   @ nxk matrix @

/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
   if ind/=0;
     al=exp(be*xl);                          @ nxk-1 matrix @
     bl=exp(h*xl);
     ah=exp(be*xh);
     bh=exp(h*xh);
     an=exp(be*xmin);
     bn=exp(h*xmin);
     ax=exp(be*xmax);
     bx=exp(h*xmax);
   endif;
   
/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
   if ind/=0;
     meanlow=al./(al+bl);                    @ nxk-1 matrix @
     meanhigh=ah./(ah+bh);
     meanmin=an./(an+bn);
     meanmax=ax./(ax+bx);
     varlow=(al.*bl)./(((al+bl)^2).*(al+bl+1));
     varhigh=(ah.*bh)./(((ah+bh)^2).*(ah+bh+1));
     varmin=(an.*bn)./(((an+bn)^2).*(an+bn+1));
     varmax=(ax.*bx)./(((ax+bx)^2).*(ax+bx+1));
   endif;
   
/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
   if ind/=0;
     fdm=meanhigh-meanlow;
     fdmmax=meanmax-meanmin;
     fdv=varhigh-varlow;
     fdvmax=varmax-varmin;
   endif;
   
   if ind/=0;
     totm2=meanc(fdm);  @ k-1X1 matrix for mean @
     totv2=meanc(fdv);  @ k-1X1 matrix for variance @
     totm2s=stdc(fdm);  @ k-1X1 matrix for mean @
     totv2s=stdc(fdv);  @ k-1X1 matrix for variance @
     totmmax=meanc(fdmmax);  @ k-1X1 matrix for mean @
     totvmax=meanc(fdvmax);  @ k-1X1 matrix for variance @
     totmmaxs=stdc(fdmmax);  @ k-1X1 matrix for mean @
     totvmaxs=stdc(fdvmax);  @ k-1X1 matrix for variance @
   endif;

   if ind/=0;
     print $ind;; totm2~totm2s~totmmax~totmmaxs;
     print $ind;; totv2~totv2s~totvmax~totvmaxs;
   endif;
  
format /rd 8,5;

alp=exp(meanx'*b1);
bet=exp(meanx'*h1);
print "mean alpha=";; alp;
print "mean beta=";; bet;
print "ave mean=";; alp/(alp+bet);
print "ave var=";; alp*bet/((alp+bet)^2*(alp+bet+1)); 

/* run ols */

{ vnam,m,bols,stb,vcols,stderr,sigma,cx,rsq,resid,dbw } = ols(0,y,x);

/* compare mse of y */
mseb=sumc((y-eyb)^2)/rows(y);
mseo=sumc((y-x*bols)^2)/rows(y);

minc(y)~maxc(y);
minc(eyb)~maxc(eyb);
minc(x*bols)~maxc(x*bols);
sumc(x*bols.>1)/rows(y);
sumc(x*bols.<0)/rows(y);

y~eyb~(x*bols);

rbro=mseb/mseo;
print "Beta is better if rb/ro<1";;rbro;

output off;

